!
!     MUNILIST   
!     Copyright (C) 1997-9  Filip Hroch, Masaryk University, Brno, CZ
!                    
! -----------------------------------------------------------------------
!

Module Mlist

contains

Subroutine munilist &
  (namfits,column,nstars,stars,outlst,min0,per,alpha,delta,long,lat, &
  sufix,status,quiet,tol)

use trajd
use traco

implicit none

character(len=*) :: namfits,sufix
character(len=1) :: outlst
integer :: column,nstars,stars(:,:),status,quiet
real(selected_real_kind(14)) :: min0,per,alpha,delta,long,lat,tol

integer, parameter :: DIM = 2
integer :: stat,BlockSize,Naxis,Naxes(DIM),PCount,GCount,BitPix
logical :: Simple,Extend

character(len=80) :: file,buf,radek

character(len=80) :: filter
integer :: i,j,y1,m1,d1,h1,min1
real(double) :: y,m,d,h,min,s,jd,ls,hcor,helcor,phase,azimut,elev,ha,airm
real :: e,x(nstars),dx(nstars),sky(nstars),dsky(nstars)
real :: is(nstars),dis(nstars)

status = -1
stat = 0
dx = 0

! Input FITS file
Call ftopen(20,namfits,0,BlockSize,stat)
Call ftghpr(20,DIM,Simple,BitPix,Naxis,Naxes,PCount,GCount,Extend,stat)
if( stat /= 0 ) then
   if( quiet == 0 ) write(*,*) "Error open file: ", namfits
   return
endif

! read date, time, exposure, filter
Call ftgkys(20,'DATE-OBS',radek,buf,stat)
Call fts2tm(radek,y1,m1,d1,h1,min1,s,stat)
y = y1
m = m1
d = d1
h = h1
min = min1
!write(*,*) y1,m1,d1,h1,min1,s

if( h1 == 0 .and. min1 == 0 ) then
   Call ftgkys(20,'UT' ,radek ,buf,stat)
   Call fts2tm(radek,y1,m1,d1,h1,min1,s,stat)
   h = h1
   min = min1
endif

!write(*,*) h1,min1,s
!read(radek,'(3(e2.0,tr1))') h,min,s

Call ftgkye(20,'EXPTIME' ,e ,buf,stat)
 if( stat /= 0 ) then
   e = 0.0
   if( quiet == 0 ) write(*,*) 'Exposure time not found. Assumed =0.'
   stat = 0
 endif

Call ftgkys(20,'FILTER' ,filter ,buf,stat)
 if( stat /= 0 ) then
   filter = '-'
   if( quiet == 0 )  write(*,*) 'Filter not found.'
   stat = 0
 endif


Call ftclos(20,stat)  
if( stat /= 0 ) return

! Evaluate Julian date at half of exposure 
jd = datjd(y,m,d + (h + (min + s/60.0)/60.0)/24.0)
jd = jd + e/2.0/86400.0

! heliocentric correction
if( alpha >= 0.0 )then
   !  approx of length of the Sun 
   ls = 279.465 + 0.985647*(jd - datjd(y,1d0,1d0))
   hcor = helcor(alpha/57.29578,delta/57.29578,ls/57.29578)
else
   hcor = 0.0
endif

! What about long filenames and extensions?
! miroslav.broz@usa.net, Dec 10th 1998

file = namfits
!i = index(namfits,'.') 
i = index(file,'fts')+index(file,'FTS')+index(file,'fits')+index(file,'FITS')
if( i == 0 ) return
!i = i + 1
if( sufix == ' ' )then
   file(i:) = 'MAT'
else
   file(i:) = sufix
endif

! Open magnitude file
open(13,file=file,status='old',iostat=j)
if( j /= 0 ) then
  if( quiet == 0 ) write(*,*) "Can't open:",trim(file)," Wrong suffix? Use 'suf='."  
  return
endif

! Read magnitude data
call rdat(13,column,nstars,stars(1:nstars,1:2),x,dx,sky,dsky)
close(13)


! Open SRT file
open(13,file=file(:i-1)//'SRT',status='old',iostat=j)
if( j /= 0 ) then
  if( quiet == 0 ) write(*,*) "Can't open:",trim(file)," Wrong suffix? Use 'suf='."  
  return
endif

! Read magnitude data
call rdat1(13,column,nstars,stars(1:nstars,1:2),real(tol),x,dx,sky,dsky)
close(13)

! Output
where( abs(dx) < 1e-4 .and. column > 3)
  x = 99.999
  dx = 9.999
endwhere

! compute intesities for stars without sky intensity
where( abs(x) < 99.0)
   is = 10**(-0.4*(x-25.0))
elsewhere
   is = -1.0
endwhere

where( abs(dx) < 9.0)
   dis = 0.921034*is*dx            ! = 0.4*ln(10)*intensity*magerr
elsewhere
   dis = -1.0
endwhere


if( outlst == 'f' ) then
   ! full output
   write(*,700) int(y),int(m),int(d),int(h),int(min),int(s), &
        int(e),filter,jd,(x(i),dx(i),i = 1,nstars)
700 format (1x,i4.4,tr1,5(i2.2,tr1),i3.1,tr1,a1,tr2,f13.5,1x, 10000f7.3)

else if( outlst == 'i' ) then
   ! intensities on output
   if( alpha >= 0.0 .and. abs(lat) <= 90.0 )then
      ha = hangle(lmst(jd,long)*15.0,alpha)
      call eq2hor(ha,delta,lat,azimut,elev)
      airm = airmass(90.0-elev)
      write(*,701) int(y),int(m),int(d),int(h),int(min),int(s), &
           int(e),filter,jd,azimut,elev,airm, &
           (sky(i),dsky(i),is(i),dis(i),i = 1,nstars)
   else
      write(*,701) int(y),int(m),int(d),int(h),int(min),int(s), &
           int(e),filter,jd,(sky(i),dsky(i),is(i),dis(i),i = 1,nstars)
   endif
701 format (1x,i4.4,tr1,5(i2.2,tr1),i3.1,tr1,a1,tr2,f13.5,1x, 1p, 1000g11.5)

else
   do i = 2,nstars
      if( x(1) > 99.9 .or. x(i) > 99.9 )then
         x(i) = 99.999
         dx(i) = 9.999
      else
         ! compute diferences
!	 write(*,*) x(1),x(i),dx(1),dx(i)
         x(i) = x(i) - x(1)
         ! their errors
         dx(i) = sqrt( dx(1)**2 + dx(i)**2 )
      endif
   enddo
   if( per < 0.0 )then
      ! Output time in Julian days.
      if( nstars == 1 )then
         write(*,710) jd + hcor,(x(i),dx(i),i=1,nstars)
      else
         write(*,710) jd + hcor,(x(i),dx(i),i=2,nstars)
      endif
   else
      ! Phase
      phase = (jd + hcor - min0)/per
      phase = phase - int(phase)
      write(*,710) phase,(x(i),dx(i),i=2,nstars)
   endif
710 format(1x,f13.5,1x,1000f7.3)
 endif

status = 0

!end

!-----------------------------------------------------------------------

contains

subroutine rdat(nunit,column,n,nlin,x,dx,sky,dsky)
      
  !  Cteni jasnosti ze souboru. Sam pozna, jestli jde 
  ! o typ AP nebo jiny a podle toho cte chyby jasnosti.

  implicit none
  integer :: nunit,column,n,nlin(:,:)
  real :: xx,x(:),dx(:),sky(:),dsky(:)

  character(len=200) :: lin 
  logical :: ap
  Integer i,j,in
      
! Rozpoznani typu souboru
  ap = .false.
  read(nunit,'(a)',err=666) lin
  if( lin(1:4) == ' NL ' )then
     read(nunit,*,err=666) i
     ap =  i == 2
  endif
      
  rewind(nunit,err=666)
  j = 0
           
! Zpracovani radku
do 
   j = j + 1
   read(nunit,'(a)',end=60,err=666) lin
   in = 0
   do i = 1,n
      if( nlin(i,1) == j .and. nlin(i,2) == 0 ) then
         in = i
      endif
   enddo
   if( in /= 0 )then
!      if( len_trim(lin) < 1 ) write(*,*) "Bad line number?"
      if( ap )then
         read(lin,800,err=666) (xx,i=1,column-1),x(in)
         read(nunit,900,err=666) sky(in),dsky(in),(xx,i=3,column-1),dx(in)
         j = j + 1
      else
         read(lin,800,err=666) (xx,i=1,column-1),x(in),dx(in)
      endif
   endif
enddo
60 continue         
   return

800 format(f6.0,14f9.3) 
900 format(f13.0,2f6.2,f8.3,11f9.3)

666 if( quiet == 0 )  write(*,*) "Bad format of file."
                             
      end subroutine

subroutine rdat1(nunit,column,n,nlin,tol,x,dx,sky,dsky)
      
  !  Cteni jasnosti ze souboru. Sam pozna, jestli jde 
  ! o typ AP nebo jiny a podle toho cte chyby jasnosti.

  implicit none
  integer :: nunit,column,n,nlin(:,:)
  real :: xx,xs,ys,x(:),dx(:),sky(:),dsky(:),tol

  character(len=200) :: lin 
  logical :: ap
  Integer i,j,in,ni,xi
      
! Rozpoznani typu souboru
  ap = .false.
  read(nunit,'(a)',err=666) lin
  if( lin(1:4) == ' NL ' )then
     read(nunit,*,err=666) i
     ap =  i == 2
  endif
      
  rewind(nunit,err=666)
  j = 0
           
! Zpracovani radku
do 
   j = j + 1
   read(nunit,'(a)',end=60,err=666) lin
   in = 0
   ni = 0
   read(lin,*,iostat=ni) xi,xs,ys
   do i = 1,n
      if( ni == 0 .and. xi > 0 .and. nlin(i,2) > 0 ) then
         if( sqrt((xs - nlin(i,1))**2 + (ys - nlin(i,2))**2) < tol ) in = i
      endif
   enddo
   if( in /= 0 )then
!      if( len_trim(lin) < 1 ) write(*,*) "Bad line number?"
      if( ap )then
         read(lin,800,err=666) (xx,i=1,column-1),x(in)
         read(nunit,900,err=666) sky(in),dsky(in),(xx,i=3,column-1),dx(in)
         j = j + 1
      else
         read(lin,800,err=666) (xx,i=1,column-1),x(in),dx(in)
      endif
   endif
enddo
60 continue         
   return

800 format(f6.0,14f9.3) 
900 format(f13.0,2f6.2,f8.3,11f9.3)

666 if( quiet == 0 )  write(*,*) "Bad format of file."
                             
 end subroutine rdat1

end Subroutine munilist

end Module Mlist
